#delimit;
clear;


use 20161021_Working_hh.dta, clear;
keep if year>=2011 & year<=2015;
drop if statefip==11;
rename statefip FIPS;
merge m:1 FIPS using "FIPS.dta", nogen;


replace Y=(Y/1000)*CPIU2016;
replace Y=Y/sqrt(numprec);

replace totalTransfers=(totalTransfers/1000)*CPIU2016;
replace totalT=totalT/sqrt(numprec);

save 20161021_Working_eh_1014.dta, replace;

clear;
save 20161021_parametersByState_EH.dta, replace emptyok;

local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;
matrix drop _all;
use 20161021_Working_eh_1014.dta, clear;
  nl (totalTransfers = {a} + {b1}*exp({b2}*Y)) if state_al=="`i'", initial(a 1 b1 2 b2 -.1) vce(gnr);
  matrix B=e(b);
  matrix V=e(V);
  matrix v=(V[1,1], V[2,1], V[2,2], V[3,1] ,V[3,2], V[3,3]);
  
	
  clear;
  svmat B;
  svmat v;
  gen state_alpha="`i'";
  append using 20161021_parametersByState_EH.dta;
  save 20161021_parametersByState_EH.dta, replace;

};


use "FIPS.dta", clear;
sort state_alpha;

merge 1:1 state_alpha using 20161021_parametersByState_EH.dta ;
drop if _merge~=3;
drop _merge;

rename B1 A;
rename B2 B1;
rename B3 B2; 
rename v1 V_A;
rename v2 C_AB1;
rename v3 V_B1;
rename v4 C_AB2;
rename v5 C_B1B2;
rename v6 V_B2; 


sort state;
save 20161021_parametersByState_EH, replace;



*>>>>A  Constructing Matrices;



use 20161021_parametersByState_EH, clear;

matrix drop _all;

mkmat A B1 B2, matrix(BETA) rown(state_alpha);



mkmat V_A C_AB1 V_B1 C_AB2 C_B1B2 V_B2  , matrix(COV) rown(state_alpha);



local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;

    matrix BETA_`i'= BETA["`i'", 1..3];
    matrix COV_`i'=COV["`i'", 1..6];


        
};


*>>>>B  Generating Simulated Parameters;

drop _all;
save "20161021_Simulations_eh.dta", replace emptyok;


local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;


foreach i of local list  {;


    drawnorm ALPHA BETA_1 BETA_2, n(1000) means(BETA_`i') cov(COV_`i') cstorage(lower) seed(00001) clear;
    generate state_alpha="`i'";
    
	append using "20161021_Simulations_eh.dta";
    save "20161021_Simulations_eh.dta", replace;

 
    
};

*>>>>C  Generating Distribution of TAU and R;


use 20161021_Simulations_eh.dta, clear;
sort state;

merge m:1 state_alpha using 20161021_parametersByState_EH; 


local list "AK AL AR AZ CA CO CT DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS
MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY" ;

gen PSI_1=.5*(19.790/sqrt(3));
gen PSI_2=1*(19.790/sqrt(3));
gen PSI_3=1.5*(19.790/sqrt(3));
gen PSI_4=2*(19.790/sqrt(3));

forvalues p=1/4{;

generate TAU_`p'=.7*PSI_`p';

foreach i of local list  {;

local j=1;

while `j'  < 100  {;
    generate TAU_`p'_`i'_`j'= TAU_`p' +`j'*.1 if (state_al=="`i'") ;
	generate TEST_`p'_`i'_`j'=(PSI_`p'-TAU_`p'_`i'_`j') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`j')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`j' if TEST_`p'_`i'_`j'<.1 & TEST_`p'_`i'_`j'>0;
	
	local j=`j'+1;
	};



local k=1;

while `k'  < 100  {;
    replace TAU_`p'_`i'_`k'= TAU_`p' - `k'*.01 if (state_al=="`i'") ;
	replace TEST_`p'_`i'_`k'=(PSI_`p'-TAU_`p'_`i'_`k') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`k')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`k' if TEST_`p'_`i'_`k'<.01 & TEST_`p'_`i'_`k'>0;
	
	local k=`k'+1;
	};


local l=1;

while `l'  < 100  {;
    replace TAU_`p'_`i'_`l'= TAU_`p' - `l'*.001 if (state_al=="`i'") ;
	replace TEST_`p'_`i'_`l'=(PSI_`p'-TAU_`p'_`i'_`l') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`l')) if (state_al=="`i'") ;
	replace TAU_`p'=TAU_`p'_`i'_`l' if TEST_`p'_`i'_`l'<.001 & TEST_`p'_`i'_`l'>0;
	
	local l=`l'+1;
	};



drop TEST* TAU_`p'_*;

};

generate R_`p'= ((A*TAU_`p' + (B1/B2)*exp(B2*TAU_`p') - (B1/B2))  + (PSI_`p'^2 - PSI_`p'^2/2) - (PSI_`p'*TAU_`p' - TAU_`p'^2/2)) /(PSI_`p'^2 - PSI_`p'^2/2);


};
                                                                                  

/* NB Corrected equation */;     


#delimit;
collapse A B1 B2  V_A V_B1 V_B2 TAU* PSI* R* C* (sd) SD_R_1=R_1 SD_R_2=R_2 SD_R_3=R_3 SD_R_4=R_4 SD_TAU_2=TAU_2, by(FIPS); 

merge 1:1  FIPS using FIPS.dta, nogen;



save "20161021_State_Parameters_w_Var_EH.dta", replace;

use "20161021_State_Parameters_w_Var_EH.dta", clear;

rename * *_EH;
rename FIPS_EH FIPS;
merge 1:1 FIPS using 20161021_State_Parameters_w_Var.dta;

twoway scatter R_2_EH R_2, m(i) mlab(state_alpha) mlabs(small) mlabp(0) 
 ytitle("Poverty Relief Ratio (R -- Equivalent Household)", size(medsmall)) 
  xtitle("Poverty Relief Ratio (R)", size(medsmall)) 
  yscale(range(0.3/0.7)) 
  xscale(range(0.3/0.7)) 
  ylabel(.3(.1).7, labsize(small)) 
  xlabel(.3(.1).7, labsize(small)) xsize(4) ysize(4)
  || function y=x, range(.3 .7) legend(off)
	;
	
graph save 20161021_Comparison_EH.gph, replace;

graph export 20161021_Comparison_EH.pdf, replace fontface("Helvetica");




